#Changed?
# Establish directories
# NB BRT Source Functions must be located in the working directory

in.dir = '/home1/99/jc152199/brt/data/'
setwd(in.dir)
out.dir = '/home1/99/jc152199/brt/data/'

# Load the gbm library to perform Boosted Regression Tree Analysis
necessary=c('gbm','SDMTools','random','sp','rgdal')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Load source code for BRT functions
source('brt.functions.R.cjsedit.r')

# Read in data to model
model.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv',header=T)

	# Create a sequence of random numbers between 1 and nrow(sub.data)
	#rn = randomNumbers(n=round(nrow(model.data)/3,0), min=1, max=nrow(model.data), col=1)
	#rn = unique(rn)
	# Create a training dataset (1/3 the size of sub.data) and testing data set (2/3 the size of sub.data)
	#train.data = model.data[c(rn),]
	#test.data = model.data[-c(rn),]

# Preferentially use gbm() before gbm.step()
	# NB - gbm.steop() doesn't work if .Random.seed function is called prior to running the model
	#
brt1 = gbm.step(data=model.data,gbm.x = c(4:8,11),gbm.y = 3,family = "gaussian",tree.complexity = 4,learning.rate = 0.01,bag.fraction = 0.5, n.folds=10)

brt = gbm(micro_max~AWAP_max+coastdist+fpcmean+fpcvar+roaddist+solar, data=model.data, distribution="gaussian", cv.folds = 10, shrinkage = 0.001)

 	for(shrinkage in c(0.001,0.0025,0.005,0.0075,0.01)) 
	
	{
 
	brt = gbm(micro_max~AWAP_max+coastdist+fpcmean+fpcvar+roaddist+solar, data=model.data, distribution="gaussian", cv.folds = 10, shrinkage = shrinkage,n.trees = 1000)
 
	out = rbind(out,data.frame(shrinkage=shrinkage,t(brt$oobag.improve)))
 
	}
	
	ylimits = range(out[,-1])
	xvals = 1:1000
	plot(c(1,1000),ylimits,type='n')
	
	for(ii in 1:nrow(out)) 
	
	{
	
	lines(x=xvals,y=out[ii,-1],col=rainbow(nrow(out))[ii])
 
	}
 
dev.off()

# gbm.perf won't test for the optimal number of iterations with cross-validation unless cv.folds is specified in the gbm call
# I think gbm.step won't work with 

check = gbm.perf(brt2, plot.it = TRUE, oobag.curve = FALSE, overlay = TRUE, method='cv')

# The partial dependency plots should be located inside the loop and written out for each model as well
# Plot fitted functions to determine importance of individual predictors

pdf(file='brtplot1')
par(mfrow=c(3,4))
gbm.plot(brt1, n.plots=6, write.title = F)
dev.off()

# Time to make spatial predictions, establish ASCII dir and read in data
# Remember that to make spatial predictions in BRT all grid data must be the same size and spatial resolution
# This means awapmax.asc will need to be resampled at 100 m resolution in UTM

ascii.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/'

fpcmean.asc = read.asc(paste(ascii.dir,'STATIC/','fpcmeanwtplusbuffer250m.asc',sep=''))
fpcvar.asc = read.asc(paste(ascii.dir,'STATIC/','fpcvarwtplusbuffer250m.asc',sep=''))
coasdtdist.asc = read.asc(paste(ascii.dir,'STATIC/','coastdistwtplusbuffer250m.asc',sep=''))
disttoroad.asc = read.asc(paste(ascii.dir,'STATIC/','disttoroadwtplusbuffer250m.asc',sep=''))
rad015.asc = read.asc(paste(ascii.dir,'SOLAR/rad015.asc',sep=''))
awapmax.asc = read.asc.gz('/home1/99/jc152199/MicroclimateStatisticalDownscale/SpatialTMax/AWAP/Tmax/2007/tmax.2007010120070101.grid.gz')

# Establish row/column positions for 250 m grid

base.asc = fpcmean.asc
base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T)) 

# Add east and north to this data frame of row/column positions

base.pos$east = getXYcoords(base.asc)$x[base.pos$row]
base.pos$north = getXYcoords(base.asc)$y[base.pos$col]

# Convert east-north to lat-long
### EPSG 32755 WGS1984 Lat Long
### EGSG 4326 WGS1984 UTM

tout = as.data.frame(spTransform(SpatialPoints(cbind(base.pos[3:4]),proj4string=CRS("+init=epsg:32755")), CRS("+init=epsg:4326")))
base.pos$lon = tout[,1]
base.pos$lat = tout[,2]

#Intersect these positions with awapmax.asc

awapmax.pos = data.frame(base.pos$col, base.pos$row, extract.data(cbind(base.pos$lon,base.pos$lat),awapmax.asc))

# Write AWAP data back into base.asc, export as an ASCII, then reimport and check

base.asc[cbind(base.pos$row,base.pos$col)]=awapmax.pos[,3]

write.asc(x=base.asc, file=paste(ascii.dir,'STATIC/AWAPmaxwtplusbuffer250m.asc',sep=""))

awapmax.asc = read.asc(paste(ascii.dir,'STATIC/AWAPmaxwtplusbuffer250m.asc',sep=''))

# Making spatial predictions from a BRT model
# Start by making a vector of the names of the ASCII grids
# Then rename the columns from model data to the same as the name of the ASCII file
# Then name the variables for the BRT model the same as both the model and ASCII grid names

grid.names = c('AWAPmaxwtplusbuffer250m.asc','coastdistwtplusbuffer250m.asc','fpcmeanwtplusbuffer250m.asc','fpcvarwtplusbuffer250m.asc','disttoroadwtplusbuffer250m.asc','rad015.asc')

names(model.data)[c(4:8,11)]=grid.names

variable.names <- c(names(model.data)[c(4:8,11)])

# Set the working directory the location where the ASCII grids are stored or the scan function won't work

setwd('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/')

for(i in 1:length(grid.names))
	
	{
	
	assign(variable.names[i],scan(grid.names[i], skip=6, na.string = "-9999"),pos=1)

	}

# Bind all the grid data into a single dataframe
	
predict.data = data.frame(AWAPmaxwtplusbuffer250m.asc,coastdistwtplusbuffer250m.asc,fpcmeanwtplusbuffer250m.asc,fpcvarwtplusbuffer250m.asc,disttoroadwtplusbuffer250m.asc,rad015.asc)
	
# Run the spatial prediction function

gbm.predict.grids(brt1, predict.data, want.grids = T, sp.name = "micro_max",pred.vec = rep(-9999,2095980), filepath = "/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/", num.col = 1086, num.row = 1930, xll = 254162.1, yll = 7812803, cell.size = 250, no.data = -9999, plot=T)

# Create an image of the ASCII

micro_max.asc = read.asc(paste(ascii.dir,'micro_max.asc',sep=''))

png(paste(in.dir,'micro_max.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(micro_max.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(heat.colors(9)[9:1]))
      
      dev.off()
	  
train.data$test = rep(1,nrow(train.data))






